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Abstract 

o 

■ We have undertaken the simulation of hydrodynamic flows with bulk Lorentz factors 
'• in the range 10^-10^. We discuss the application of an existing relativistic, hydro- 
I dynamic primitive- variable recovery algorithm to a study of pulsar winds, and, in 

■ particular, the refinement made to admit such ultra-relativistic flows. We show that 
an iterative quartic root finder breaks down for Lorentz factors above 10^ and em- 

O I ploy an analytic root finder as a solution. We find that the former, which is known 

' to be robust for Lorentz factors up to at least 50, offers a 24% speed advantage. We 

Ci ■ demonstrate the existence of a simple diagnostic allowing for a hybrid primitives 

recovery algorithm that includes an automatic, real-time toggle between the itera- 
tive and analytical methods. We further determine the accuracy of the iterative and 
hybrid algorithms for a comprehensive selection of input parameters and demon- 
5^ I strate the latter's capability to elucidate the internal structure of ultra-relativistic 

plasmas. In particular, we discuss simulations showing that the interaction of a 
light, ultra-relativistic pulsar wind with a slow, dense ambient medium can give rise 
to asymmetry reminiscent of the Guitar nebula leading to the formation of a rela- 
tivistic backflow harboring a series of internal Shockwaves. The Shockwaves provide 
thermalized energy that is available for the continued inflation of the PWN bubble. 
In turn, the bubble enhances the asymmetry, thereby providing positive feedback 
to the backflow. 
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1 Introduction 



Hydrodynamic simulations have been widely used to model a broad range 
of physical systems. When the velocities involved are a small fraction of the 
speed of light and gravity is weak, the classical Newtonian approximation 
to the equations of motion may be used. However, these two conditions are 
violated for a host of inte r esting scenarios, including, for exa mple, heavy ion 



collis ion systems (IHirand . 120041 ). r elativis t ic las er systems (jPelettrez et al. 



20051 ) , and many from astrophysics (jibanezl . l2003l . and references therein) , that 
call for a fully relativistic, hydrodynamic (RHD) treatment. The methods of 
solution of classical hydrodynamic problems have been successfully adapted 
to those of a RHD nature, albeit giving rise to significant complication; in 
particular, the physical quantities of a hydrodynamic flow (the rest-frame 
mass density, ra, pressure, p, and velocity, v) are coupled to the conserved 
quantities (the laboratory-frame mass density, R, momentum density, M, and 
energy density, E) via the Lorentz transformation. The fact that modern RHD 
codes typically evolve the conserved quantities necessitates the recovery of the 
physical quantities (often referred to as the "primitive variables") from the 
conserved quantities in order to obtain the flow velocity. Thus, the calculation 
of the primitives from the conserved variables has become a critical element 
of modern RHD codes (IMartf and Miilleii . l2003l ). Indeed, this is an active area 
of research with significant attention given t o the case of gener al relativistic, 
magnetohydrodynamic (GRMHD) case (e.g. iNoble et all. l2006h and varying 



equat ions of state wi thin the context o f RMHD (e.g. iMignone and McKinneyi . 
20071 ) and RHD (e.g. iRyu et al.l . l2006l ). This work is concerned with the RHD 
case for a fixed adiabatic index and so we refer the reader to the above- 
mentioned papers for a discussion of those studies. 



In this paper, we present a method for recovering the primitive variables 
from the conserved quantities representing special relativistic, hydrodynamic 
(SRHD) flows with bulk Lorentz factors (7 = (1 — f^)~^/^, where v is the bulk 
flow velocity - the speed of light is normalized to unity throughout this paper) 
up to 10^. We started with a module f rom an existing SRHD code used to sim- 
ulate flows with 7 <50 as described in lDuncan and Hughes! (Il994j ). Admitting 
flows with such ultra-relativistic Lorentz factors as 10^ required significant re- 
finement to the method used in the existing code to calculate the flow velocity 
from the conserved quantities. In particular, such extreme Lorentz factors lead 
to severe numerical problems such as effectively dividing by zero and subtrac- 
tive cancellation. In §2, we discuss the formalism of recovering the primitives 
within the context of the Euler equations. In §3, we elucidate the details of 
the refinement to this formalism necessitated by ultra-relativistic flows. We 
present the refined primitives algorithm in §4 and our application in §5. We 



(P. A. Hughes). 
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discuss our results and conclusions in §6 and §7, respectively. 



2 Recovering the primitive variables from R, M, and E 



In general, recovering the primitives from the conserved quantities reduces to 
solving a quartic equation, Q{v) = 0, for the flow velocity in terms of R, M, 
and E. Implementation typically involves a numerical root flnder to recover 
the velocity via Newton- Raphson iteration which is very eflicient and provides 
robustness because it is straightforward to ensure that the computed velocity 
is always less than the speed of light. This is a powerful method that is in- 
dependent of dimensionality and symmetry. The latter point follows directly 
from the fact that symmetry is manifest only as a source term in the Euler 
equations and does not enter into the derivation o{Q{v) (see the axisymmetric 
example below). Dimensional generality arises because regardless of the coor- 
dinate system. 



one may always write M = M^. , where the M^^ are the 
components of the momentum-density vector along the orthogonal coordinates 
Xj. In the case of magnetohydrodynamic (MHD) flows, there are, of course, ad- 
ditional considerations. However, non-magnetic (RHD) simula tions sti l l have 
a signiflcant role to play ir i astrophysics, e.g. extragala ctic jets (jHughesl . l2005l ) 



and pulsar wind nebulae ( van der Swaluw et al. . 20041 ) 



As an example, consider the case of the axisymmetric, relativistic Euler equa- 
tions, which we apply to pulsar winds. This type of formalism enjoys diverse 
application, in both specia l and general r elatiy istic settings, from 3D simula- 



tions of extragalactic jets (IHughes et al 



of gamma-ray bursts (jZhang et al.l. l2003r) an d the collapse of massive stars to 



20021 ). to theories of the generation 



neutron stars and black holes (jShibatal . l2003l ). In cylindrical coordinates p and 
z, and deflning the evolved- variable, flux, and source vectors 



U = {R,M„M,,Ef, 
FP = {RyP, My + p, My, {E + p^f, 

F' = {Rv', My, My + p,{E + pyf, 

S={0,p/p,0,Of, (1) 
the Euler equations may be written in almost-conservative form as: 



The pressure is given by the ideal gas equation of state p = (F — l)(e — n), 
where e and F are the rest-frame total energy density and the adiabatic index. 
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Note that the velocity and pressure appear exphcitly in the relativistic Euler 
equations, in addition to the evolved variables, and pressure and rest density 
are needed for the computation of the wave speeds that form the basis o f 
typical numerical hydrodynamic solvers, such as that due to iGodunovl (Il959l ). 
We obtain these values by performing a Lorentz transformation where the 
rest-frame values are required: 



R = 'yn, 



"1/2 



(2) 



where = {v^Y + {v^f and = -i^{e + pY[{vPf + [v^'f] = j^^e + p)'^v'^. 
When the adiabatic index is constant, combining the above equations with the 
equation of state creates a closed system which yields the following quartic 
equation for v in terms of Y = M/E and Z = R/E: 



Q{v) = (r - iy{Y' + z'y - 2r(r - i)Yv 



+ 



r' + 2(r- (r- 1)2^2 



- 2rYv + Y^ 



0. 



(3) 



Component velocities, and the rest-frame total energy and mass densities are 
then given by: 



M 

e = E - Mpvp - M,v\ 



n - 



R 



3 Refinement of the root finder to admit ultra-relativistic flows 



A particular implementation of the a bove has been previously a pplied to rel- 



ativistic galactic jets with 7 < 50 (IDuncan and Hughed . [l99J). The ultra- 
relativistic nature of pulsar winds necessitated an investigation of the behav- 
ior of the primitives algorithm upon taking 7 >> 1. We found that, beyond 
7 ~ 10^, the algorithm suffers a severe degradation in accuracy that worsens 
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with increasing Lorentz factor until complete breakdown occurs due to the 
failure of the Newton-Raphson iteration process used to calculate the flow 
velocity. 

The problem lies in the shape of the quartic, Q(f ), one must solve to calculate 
the primitive variables. The quartic equation as derived using the velocity as 
a variable exhibits two roots for typical physical parameters of the flow (see 
Fig. 1). In general, for 7 < 10^, the two roots are sufficiently separated on the 
velocity axis such that the Newton-Raphson (N-R) iteration method converges 
to the correct zero very quickly and accurately (for Y < 0.9 and Z > 10~^, 
corresponding to 7 < 2, the roots approach each other sufficiently such that 
the incorrect root is selected; see §4.3). In fact, N-R iteration can be so efficient 
that it is more desirable to use this method than it is to calculate the roots 
of the quartic analytically (see §4.2). However, as the Lorentz factor of the 
flow increases, the roots move progressively closer together and the minimum 
in Q{v) approaches zero. Eventually, the minimum equals zero to machine 
accuracy which causes dQ/dv = to machine accuracy resulting in a divide 
by zero and the Newton-Raphson method fails (see Fig. 2). 

A simple and highly effective solution (see §4.3 for details) is to rewrite the 
velocity quartic, Q{v) (Eqn. 3), in terms of the Lorentz factor (i.e. make the 
substitution v"^ = 1—7^^) to obtain the quartic equation in 7 (recall Y = M/E 
and Z = R/E): 



Q(7) = T\i - r^)^ _ 2r(r - 1)^7^ 
+ [2r(r - 1)^2 + (F - ifz"^ - F^" 

+ 2F(F - 1)^7 - (F - l)2(y2 + 



7 



0. 



(4) 



As Fig. 1 exemplifies, ^(7) exhibits a single root for the physical range 7 > 
1. However, Newton-Raphson iteration also fails in this case at high Lorentz 
factors because of the steepness of the rise in Qipf) through the root. Thus, we 
are forced to use an analytical method of solving a quartic. Below, we discuss 
our implementation. 



3.1 Solving a quartic equation 



We use the prescription due to iBronshtein and SemendyayevI (119971 ) in order 
to analytically solve for the roots of a quartic. We chose this method because 
it provides equations for the roots of the quartic that are the most amenable 
(of the methods surveyed) to integration into a computational environment. 
In order to provide a complete picture of our m ethod, which includes steps not 
found in IBronshtein and SemendyayevI (1l997l ). we reproduce some sections of 



5 



Fig. 1. The left-hand plots show the shape of the Lorentz factor quartic over a run 

of Lorentz factors for a mildly rclativistic flow (70 = 1.5) and an ultra-relativistic 
flow (70 = 10^). The right-side plots show the shape of the velocity quartic over a 
run of velocity for a mildly relativistic flow (70 ~ 1.5) and a highly (but not ultra) 
relativistic flow (70 ~ 10^). The crosses mark the location of the physical root. Prom 
the plot in the lower right, one can see the onset of the zero derivative problem as 
the roots are not distinguishable from each other or the local minimum even on a 
scale of 10~^^, which begins to encroach on the limit of 8-byte accuracy. 
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that text. We precede as follows. 
Given a quartic equation in x: 

a^x'^ -\- a'iX^ -\- a2X^ -\- aix -|- oq = 0, a„ G 9?, 04 7^ 0, (5) 

normalizing the equation (dividing by 04) and making the substitution y ~ 
X + results in the reduced form: 

y^ + Py^ + Qy + R^Q, 
where, defining a„ = 0^/04: 

3 2 
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Table 1 

The dependence of the sohitions to the parent qnartic on the sohitions to the cubic 
resolvent. 



Solutions of the cubic resolvent Solutions of the quartic equation 

all real and positive all real 

all real, one positive two complex conjugate (cc) pairs 

one real, one cc pair two real, one cc pair 



These coefficients allow the definition of the cubic resolvent: 

+ 2Pu^ + _ /^R)u - = 0, (6) 



upon whose solutions the solutions of the original quartic (Eqn. 5) depend. 
The product of the solutions of the cubic resolvent is U1U2U3 — Q"^ (Vieta's 
theorem), which clearly must be positive. The characteristics of the quartic's 
roots depend on the nature of the roots of the cubic resolvent (see Tab. 1). 

Given the solutions of the cubic resolvent Ui, U2, and U3, the solutions of the 
quartic (Eqn. 5) are 
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(7) 
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3.2 Solving a cubic equation 



The equations of the previous section reduce the problem of solving a quartic 
equation to that of solving a cubic equation (i.e. the cubic resolvent of Eqn. 6). 



Once again following iBronshtein and SemendyayevI (119971 ) (note the similarity 
to the method in the previous section), given a cubic equation: 



hu^ + h2U^ + hiu + 6o = 0, 6„ e 3?, 63 ^ 0, 



normalizing the equation and making the substitution v = u + &2/3&3 results 
in the reduced form: 

+ + g = 0, 
where, defining 6„ = hn/h^'- 



These coefficients allow the definition of the discriminant. 




upon which the characteristics of the solutions of the cubic equation depend 
(see Tab. 2). 

Given p, g, and Cardando's formula for the reduced form of the cubic leads 
to the solutions of the original cubic (Eqn. 8): 



Ui = 


s + t 


b2 

363' 


U2 = 




+ t)- 


U3 = 




+ t)- 


where: 







b2 . .^3. ,^ 

— — h I s — t), 

363 2 ^ ^' 

(10) 
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Table 2 

The dependence of the solutions of a cubic equation on the sign of the discriminant 
(assuming a real variable). 



D 



Solutions of the cubic equation 



positive one real, one complex conjugate pair 



negative 



all real and distinct 



= 



all real, two (one, if p = g = 0) distinct 



s 



t 



li D < 0, the the cubic has three real roots, subject to the following two 
subcases, and the four real roots of the quartic follow directly from Eqn. 7. If 
D — 0, then s = t and the cubic has three real solutions that follow directly 
from Eqn. 10 from which one can see that two are degenerate. If D < 0, 
the cubic has three distinct real roots. Obtaining these sohitions via Eqn. 10 
requires intermediate complex arithmetic. However, this may be circumvented 
by making the substitutions: 




COS(0) 



2r' 



in which case the solutions of the cubic (Eqn. 8) are: 




(11) 
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If D > 0, then the cubic has one real root and a pair of complex conjugate 
roots and the quartic has two real roots and a pair of complex conjugate roots 
(see Tab. 1). Finding the roots of the quartic involves intermediate complex 
arithmetic which may be circumvented as follows. Defining: 



Eqn. 10 may be rewritten as: 



Ui=S + t- 

303 

U2 = R + iC, 
U3 = R — iC. 

Next, we have 1*2,3 — VR^ + (7^6=*=*^/^. We then obtain the roots of the quartic 
from Eqn. 7: 



Note that xi and X2 are the two real solutions. 



4 The refined primitives algorithm 



Using the method above we created a SRHD primitive algorithm called "REST_ 
FRAME". Given the speed advantage of the iterative root finder (see §4.2), 
it a desirable choice over the analytical method within its regime of ap- 
plicability, i.e. for low Lorentz factors. As Fig. 2 shows, the iterative root 
finder is accurate to order 10~^ (see §4.3) for a sizable region of parameter 
space including all R/E above the diagonal line between the points (0, —7) 
& (9, 0) in the log(i?/^) vs. - log(l - M/E) plane (i.e. for log(i?/^) > 
— (7/9) X log(l — M/E) — 7). Therefore, for a given M/E and R/E, we check 
if this inequality is true; if (not) so, we call the (analytical) iterative root finder 
(see §4.1). 
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4-1 Pseudo-code 



REST_FRAME calculates the primitive variables given the conservative vari- 
ables and the adiabatic index as represented in the following pseudo-code (note 
this is a 2D example): 

PROCEDURE REST_FRAME 

RECEIVED FROM PARENT PROGRAM: Y, Z 

RETURNED TO PARENT PROGRAM: 7, v, C 
Comment: recall Y = M/E and Z = R/E 
Comment: C is returned < for code failures 

GLOBAL VARIABLE: P 

SET VALUE OF rriunderfiow 

SET VALUE OF vtoi 
Comment: determines iterative method velocity accuracy 
Comment: we set Vtoi = 10-^ 10^^^ lO^^^^ 10^^^ 

Comment: for — log(l — Y)< 8.3, < 10.3, < 12.3, otherwise, respectively 
SET M = ^MJ + M2 

IF M < munderflow THEN 
f = 0, 7 = 1 

Comment: avoids code failure if v is numerically zero 
ELSE 

TEST FOR UNPHYSICAL PARAMETERS 
IF PASSED, SET C NEGATIVE AND RETURN 
IF log(Z) > -(7/9) X log(l - y) - 7, THEN 
Comment: check to see if input parameters are within the acceptable 
Comment: accuracy region of the iterative routine 

CALL ITERATIVE_QUARTIC(y, Z, vtou v,C) 
Comment: updates Vn-i to Vn using n cycles of Newton-Raphson iteration 
Comment: returns v = f„ when \vn — Vn-i\ < Vtoi 

IF C < 0, THEN 
Comment: this means the iteration failed to converge 
RETURN 
ELSE 

END IF 
ELSE 

CALL ANALYTICAL_QUART1C(F, Z, 7) 
Comment: calculat es 7 usi ng analytical solution - see below 

END IF 

END IF 

END PROCEDURE REST_FRAME 



11 



PROCEDURE ANALYTICAL_QUARTIC 
Comment: see §3.1 for equations 

RECEIVED FROM PARENT PROGRAM: F, Z 
RETURNED TO PARENT PROGRAM: 7 
GLOBAL VARIABLE: F 
as = 2F(F-l)Z(F-2 + l) 

02 = (F2 - 2F(F - 1)^2 _ (r - l)2z2)(r-2 + 1) 
a\ = —as 

ao = {T-iy{Y^ + Z^){Y-^ + l) 

54 = 1 + y2 _ _ Q2 

Comment: coefficients recast to counter subtractive cancellation - see §4.3 

NORMALIZE COEFFICIENTS TO 04 
Comment: e.g., a^N = a^/a^ 

CALCULATE CUBIC RESOLVENT COEFFICIENTS 

CALCULATE DISCRIMINANT, D 

IF D <0 THEN 

WRITE ERROR MESSAGE AND STOP 
Comment: exploration suggests D < is unphysical but formal proof is elusive 
Comment: thus, we leave D <0 uncoded with a error flag just in case 

ELSE 

Comment: D > ^ Q{l) has 2 real roots (see Tab. 1 & 2) 
CALCULATE ROOTS OF CUBIC RESOLVENT 

Comment: the cubic has one real root and a pair of complex conjugate roots 
IF REAL ROOT < 0, SET REAL ROOT = 

Comment: the real root cannot be less than zero analytically 

Comment: numerically, however, it can have a very small negative value 
CALCULATE THE TWO REAL ROOTS OF THE QUARTIC 
TEST FOR TWO OR NO PHYSICAL ROOTS 
IF PASSED, WRITE ERROR MESSAGE, AND RETURN 
IF FAILED, SET 7 = PHYSICAL ROOT 
END IF 

END PROCEDURE ANALYTICAL_QUARTIC 



4.2 Code timing 

Using the Intel Fortran library function CPU_TIME, we calculated the CPU 
time required to execute 5x10^ calls to REST_FRAME for Y = 0.9975 & Z = 
lxlO~^ (7 ~10) using the Newton- Raphson iterative method with Q{v) and 
8-byte arithmetic, and the analytical method with Q(7) and both 8-byte & 
16-byte arithmetic (we investigated the use of 16-byte arithmetic due to an 
issue with subtractive cancellation - see §4.3). The CPU time for each of these 
scenarios was 29.5, 36.5 (averaged over ten runs and rounded to the nearest 
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half second) , and ~11650 seconds (one run only), respectively. This indicates 
that while using the 8-byte analytical method is satisfactory, it is advantageous 
to use the iterative method when Lorentz factors are sufficiently low, and that 
the use of 16-byte arithmetic is a nonviable option. This result is not surprising 
as the accuracy of Newton- Raph son iteration improve s by a pproximately one 
decimal place per iterative step ( Duncan and Hughed . 1994) and the relative 
inefficiency of 16-byte arithmetic is a known issue (e.g. 



Perret-Gallixl . 120061 ) 



4-3 Solver Accuracy 



The input parameters for our primitives algorithm are the ratios of the laboratory- 
frame momentum and mass densities to the laboratory-frame energy density 
(recall Y = M/E and Z = R/E) both of which must be less than unity in or- 
der for solutions of Eqn. 2 to exist. In addition, the condition Y'^ + Z'^ < 1 must 
be met. Along with the fact that Y and Z must also be positive, this defines 
the comprehensive and physical input parameter space to be < F, Z < 1 
such that Y"^ + Z"^ < 1 (we identify a particular region of parameter space 
applicable to pulsar winds in the next section) . We tested the accuracy of our 
iterative and hybrid primitives algorithms within this space as follows. 

First, as we are most interested in light, highly relativistic flows (i.e. Z small 
and Y close to unity), to deflne the accuracy-search space we elected to use the 
quantities — log(l — Y), which for values greater than unity gives 0.9 < F < 1, 
and log(Z), which for values less than negative unity gives Z <ti 1. We selected 

< — log(l — y) < 13 and —13 < log(Z) < corresponding to Lorentz factors 
(7) between 1 and 2 x 10^. We chose a range with a maximal 7 slightly above 

1 X 10^ in order to completely bound the pulsar wind nebula parameter space 
deflned in the next section. 



Choosing a relativistic equation of state r= 4/3 and using 1300 points for both 
— log(l — Y) and log(Z), we tested the accuracy of REST_FRAME by passing 
it Y and Z, choosing E = 1, and using the returned primitive quantities to 
derive the calculated energy density Ec, and calculating the difference |1 — 
Ec/E\ = 6E/E. We chose this estimate of the error because 6E/E ^ S'j/'-f 
and 57/7 is tied to the accuracy of the numerical, hydrodynamic technique 
(see the flnal paragraph in this section). 

Our results for the Newton-Raphson (N-R) and hybrid methods are given in 
Figs. 2 & 3 which show where the accuracy is of order at least 10~^, at least 
10~^, worse than 10~^, failure, and unphysical input {Z"^ > 1 — Y"^), respec- 
tively. We chose an accuracy of order 10~^ as the upper cutoff because N-R 
iteration returns accuracies on this or der for 7 < 50 and re l ativis tic, hydro- 



dynamic simulations of galactic jets by lDuncan and Hughes! (119941 ) produced 
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robust results for Lorcntz factors of at least 50 using N-R iteration. An addi- 
tional result of interest is that the ultra-relativistic approximation for v (i.e. 
taking = thereby reducing Q{v) = to a quadratic equation) manages an 
accuracy of at least 10~^ for a large portion of the physical Y — Z plane (see 
Fig. 4). 

Fig. 2 shows the accuracy of the N-R iterative method. There are several 
noteworthy features. First is the presence of a sizable region corresponding 
to 7 < 500 within which accuracy is generally significantly better than 10~^. 
Second is that N-R iteration is unreliable due to sporadic failures for increasing 
Lorcntz factors until accuracy becomes unacceptable or the code fails outright 
due to divide by zero (see §3) or non-convergence within a reasonable number 
of iterations. In addition, though N-R iteration has been widely established 
as the primitives recovery method of choice for fiows with Lorentz factors less 
than order 10^, we found that for a subset of parameters, corresponding to 7 < 
2, our N-R algorithm suffered an unacceptable degradation in accuracy. The 
key to this problem lies in the how the flow velocity {v) is initially estimated 
for the first iterative cycle as follows: 

Fig. 2. The accuracy (estimated as SE/E) of the Newton- Raphson (N-R) iterative 
primitives algorithm where white, light grey, medium grey, dark grey, and hatched 
regions correspond, respectively, to an accuracy of order at least 10~^, at least 
10~^, worse than 10~^, failure, and unphysical input {E? /E"^ > 1 — M"^ /E"^). Note 
that the Lorentz factor varies from order 1 at the far left to order 10^ at the far 
right. There is a sizable white region representing M/E < 0.999999 (7 < 500) and 
R/E > 5 X 10~® within which accuracy is generally significantly better than 10~^. 
N-R iteration is unreliable due to sporadic failures for all M/E and R/E such that 
R/E < 5 X 10~^ and for an ever increasing fraction of R/E > 5 x 10~® as M/E 
increases until accuracy becomes unacceptable or the code fails outright for M/E 
and R/E such that M/E > 0.999999. Faihucs are due to divide by zero (see §3) or 
nonconvergence within a reasonable number of iterations. 




2 4 6 8 10 12 

-log(1 -M/E) 
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(1) The established approach (iDuncan and Hughesl . 1 1994 ISchneider et aL . 
19931 ) is to bracket v with 



: mm(l, F + 5), 

r-^r2-4(r-i)r2 
2r(r- 1) 



(13) 



where 5 ~ 10 ^ and f^m is derived by taking the ultra-relativistic hmit 
(i.e. i? = 0) 

(2) The initial velocity is then Vi = {vmin + Vmax)/'^ + ^, where i] = {1 — 
Z){vmin - Vmax) for v-max > £ and 7] = othcrwisc (e order 10"^) 

(3) This method fails due to selection of the incorrect root when the roots 
converge. 

(4) Thus, we make a simpler initial estimate of Vi = Vmax, which guarantees 
that Vi is "uphill" from v for all physical Y — Z space and that N-R 
iteration converges on v. 

Fig. 3 shows that our hybrid algorithm REST_FRAME is accurate to at least 
10~^ for all but a smattering of the highest Lorentz factors. In fact, it is 
significantly more accurate over the majority of the physical portion of the 
Y — Z plane. The space between the parallel lines represents the PWN input 
parameters discussed in the next section. We find that multiplying (5(7) by 
(y2 _y-2-j g^j^^ rewriting the new (04) in terms of the new 02 (02) and new 

Fig. 3. The accuracy (estimated as 6E/E) of the hybrid primitives algorithm where 
white, light grey, and hatched regions correspond, respectively, to an accuracy of 
order at least 10~^, at least 10~^, and nonphysical input {R^ /E"^ > 1 — M^/E'^). 
Note that the Lorentz factor varies from order 1 at the far left to order 10^ at the far 
right. The space between the parallel lines represents PWNe input parameter space. 
The accuracy degradation at the extreme right is due to subtractive cancellation in 
the 4*'^-order coefficient of the Lorentz-factor quartic as M/E — >1. 




-log(1 -M/E) 
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oo (cio), e.g. 04 = 1 + — ao — 02, improves the accuracy somewhat, but does 
not entirely mitigate the problem. The issue of accuracy loss at large Lor entz 
factors in 8-byte primitives algorithms is a known issue (e.g. iNobld . l2003l ) for 
which we know of no complete 8-byte solution. Employing 16-byte arithmetic 
provides spectacular accuracy, but introduces an unacceptable increase in run 
time (see §4.2). 



The issue of what constitutes an acceptable error in the calculated Lorentz 
factor is decided by the fact that a fractional error in 7 translates to the same 
fractional error in p and n which are needed to calculate the wave speeds 
that forr n the bas i s of t he numerical, hydrodynamic technique, a Godunov 
scheme (iGodunovl . Il959l ) which approximates the solution to the local Rie- 
mann problem by employing an estimate of the wave speeds. We do not know 
a priori how accurate this estimate needs to be, and so procede with 8-byte 
simulations of pu lsar winds with the expectation of using shock-tube tests 
(IThompsonl . 119861 ) to validate the accuracy of the computation of well-defined 
flow structures as we approach the highest Lorentz factors. It is also notewor- 
thy that while 7 = 10^ is the canonical bulk Lorentz factor for pulsar winds, 
7 = 10^ and 10^ are still in the ultra-relativistic regime, and it may very well 
prove to be that these Lorentz factors are high enough to elucidate the general 
ultra-relativistic, hydrodynamic features of such a system. The hybrid algo- 
rithm achieves accuracies of at least 10~^ for 7 ~ 10^, which is safely in the 



Fig. 4. The accuracy (estimated as 6E/E) of the ultra-relativistic approximation 
of the flow velocity where white, light grey, medium grey, and hatched regions 
correspond to an accuracy of order at least 10~^, at least 10~^, worse than 10~^, 
and unphysical input {B? /E"^ > 1 — M^/£J^), respectively. Note that the Lorentz 
factor varies from order 1 at the far left to order 10^ at the far right. The accuracy 
degradation at the extreme right is due to the fact that the fractional error in the 
Lorentz factor is proportional to the fractional error in the velocity divided by l — v"^ 
which diverges as v — >1. 
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acceptable accuracy regime. 



5 Application to Bow-shock Pulsar Wind Nebulae 



At the end of a massive star's life, the collapse of its core to a compact object, 
i.e. a neutron star or black hole, drives a shockwave into its outer layers, 
thereby heating and ejecting them into the interstellar medium (ISM) in a 
supernova (SN) explosion. Subsequently, the shockwave overtakes the ejecta, 
expanding into the ISM, and forming a supernova remnant (SNR). Typically, a 
SN releases ~10^^ erg of mechanical energy that drives expansion of the SNR, 
sweeping up ISM material, heating it to X-ray temperatures and infusing it 
with fusion products beyond lithium. 

In a subclass of SNRs, for progenitor masses between 10 and 25 solar masses 



(e.g. iHeger et al.l . l2003l ). the compact object formed in the SN explosion is a 
rapidly-spinning, highly-magnetized neutron star surrounded by a magneto- 
sphere of charged particles. The combination of the rotation and the magnetic 
field gives rise to extremely powerful electric fields that accelerate charged 
particles to high velocities. The magnetic field interacts with the charged 
particles resulting in the spin-down of the neutron star, and the release of 
spin-down energy. A relatively small fraction of this energy is converted into 
beamed emission, manifest as an apparent pulse if the neutron star's rotation 
sweeps the beam across the Earth, leading to the designation "pulsar " . The 



bulk of the spin-down energy is converted into a pulsar wind (IMichell . Il969l ) 



which is terminated at a strong shock, downstream of which the flow is indis- 



tingu ishable from being spherically symmetric (e.g. IChatterjee and Cordes 



20021 ). The wind particles interact with the magnetic field causing them to 
emit synchrotron radiation, forming a pulsar wind nebula (PWN). The Crab 
Nebula, formed in the SN explosion of 1054 CE, is the canonical object of this 
type. The Crab exhibits pulsations from the radio, all the way up to X-rays, 
and is a prodigious source of 7-rays. 

The wind in the immediate vicinity of the pulsar is a diffuse, relativistic gas 
unlikely to be directly observable. However, the c lassic structure of f orward 



and reverse shocks separated by a contact surface (I Weaver et al.l . 119771 ) arises 
from the wind interaction with the SNR or ISM. A probe of this interaction 
is provided by optical emission from the swept-up ambient ISM, thermal X- 
ray emission from the SNR and/or the shocked ISM, and X-ray synchrotron 
emission from the s hocked wind. Furthermore, the high space velocity that 



is typical of pulsars (j Cordes and ChernofB . Il998l ) implies an asymmetric ram 



pressure on the pulsar wind from the denser ambient medium. The details of 
the morphology and of the distribution of the density, pressure, and velocity 
within the PWN depends upon the density, speed, momentum, and energy flux 
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of the pulsar wind. Thus, comparison of PWN simulations with observational 
data can provide an unparalleled method for investigating pulsar winds and, 
hence, how the surrounding medium taps the rotational energy of the pulsar. 



Pacini and Salvatil (119731 ) and lRees and GunnI (119741) pioneered t he basi c model 
of PWNe; a model furth e r deve l oped by iKennel and Coronitil (ll984alJbl') and 
Emmering and Chevalier (1987). Gaensler and Sland ( 2006 ) and Bucciantini 



( 120081 ) are excellent reviews on observational and theoretical studies of PWNe, 
respectively. For a number of reasons, a detailed, quantitative study of PWNe 
is now particularly timely. First, there is a cornucopia of high-quality data 
from space-born observatories such as the Chandra X-ray Observatory and 
XMM-Newton. Second, the total energy radiated by PWNe accounts for only 
a small fraction of the spin-down energy, leaving a large energy reservoir avail- 
able for interaction with the SNR and acceleration of ions, the partitioning of 
which is not well understood. 



Effor t s to model PWNe spa n three de cades (with seminal paperslRees and Gunnl . 
1974; Kennel and Coroniti . 1984a| jbl: lEmmering and Chevalier . 1987). While 
the case for a non-isotropic pulsar wind energy flux has long been made 
( iMichell . Il973l ). it has only been recently that a theoretical explanation of the 
mechanism behind the jet/torus structure interior to the t ermin ation shock 
has been p ut forward, and th at the predicti ons of iMichell (119731) have been 
confirmed ( Bucciantini . 20081 ). In particular, Bucciantini (2008) highlighted 
that a detailed description has been made possible by the increase in the effi- 
ciency and robustne s s of re l ativistic, numerical MH D co des stemming from th e 
work of iKomissarovl (119991). iDel Zanna et al.l (120031 ) and lGammie et al.l (120031 ). 
Simulations by iDel Zanna et al.l (120041 ) indicate that where jet formation in 
PWNe takes place is tied to where the magnetic field attains equipartition, at 
which point the magnetic filed can no longer be compressed. If this happens 
close to the termination shock, then, due to the mildly relativistic nature of 
the post-shock fiow, hoop stresses can become efficient and most of the fiow is 
diverted back toward the axis and collimated. The magnitude of the magneti- 
zation is key: if it is too small, the equipartition is reached outside the nebula, 
hoop stresses remain inefficient, and no collimation is produced. 

Other modeling, for example, that presented herein, is concerned with the 
global structure of PWNe. The enormous acceleration of the wind at the 
termination shock smears out asymmetries leaving an essentially spherically 
symmetric shocked flow to pr oduce large sc a le PW N features. In particular, 
Bucciantini et al.l (120051 ) and IVigelius et al.l (120071 ) are two recent examples 
of simul ations addressing t he str uctures that arise in bow-shock PWNe (see 



below). 'Bucciantini et al.l (120051) w e re the flrst to apply a fully-relat ivistic 



MHD code (IDel Zanna et al.l . l2003l : iDel Zanna and Bucciantinil . |2002|), and 



for an axisymmetric ge ometry, obta i ned a relativistic backflow behind the 
pulsar, as predicted by I Wang et al.l (119931 ) for PSR1929+10. However, the 
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wind Lorentz factor and pulsar velocity were 10 and 9000 km s~^, respec- 
tively, which are far from the typical values of 10^ and 500 km s~^ (indeed, 
the Guitar pulsar, the fastest known, has a transverse velocity of ~1700 km 
s~^). In addition, the paper d oes not address the "bubble" in the Guitar (see 
Fig. 5). IVigelius et al.l (120071 ) performed non- relativistic, hydrodynamic sim- 
ulations with a relax ation to cylindrical symmetry. The full 3-D FLASH code 
( iFryxell et al.l . |2000| ) was employed and an anisotropic pulsar wind, cooling 
of the shocked ISM, ISM density gradients, and ISM walls were considered. 
While the authors employed a realistic pulsar velocity of 400 km s~^, the non- 
relativistic nature of the simulations limited the Lorentz factor to order unity. 
In this section, we present fuUy-relativistic, axisymmetric, hydrodynamic sim- 
ulations of bow-shock PWNe for a realistic pulsar velocity and wind Lorentz 
factor. In particular, we address the origin of relativistic backflows leading to 
a persistent nebular bubble. 



5. 1 Bow-shock formation 



The evolution of PWNe can be broken into four broad phases: 1) free-expansion, 
2) SNR reverse shock interaction, 3) expansion in side a Sedov SNR, and 4) 



bow shock formation (for a detailed discussion, see iGaensler and Sland . 12006 



and references therein). In this work, we investigate the last stage of evo- 
lution. The time it takes fo r the pulsar to cross the SNR was obtained by 
van der Swaluw et al.l (2003): 



tcr = 1.4 X 10" 



Eq 
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(14) 



where Vpgr is the velocity of the pulsar. Once the PWN-SNR system has 
evolved to the Sedov- Taylor stage, the time elapsed is sufficiently large that 
is possib le for the pulsar to have re ached the edge of the nebula, or even 



beyond (Ivan der Swaluw et al.l . l200ll ). Thus, the pulsar escapes its original 
wind bubble, leaving behind a "relic" PWN, and traverses the SNR while 
inflating a new PWN. As the pulsar moves away from the center of th e 
remnant, the sound speed d ecreases. Following Ivan der Swaluw et al.l (119981 ) . 
van der Swaluw et al.l (120041 ) calculated the Mach number of the pulsar, M.psr, 



and found that M.psT exceeds unity after a time t = 0.5tcr, at which point 
the pulsar has travelled a distance Rpsr — 0.677Rsnr, and the nebula is de- 
formed into a bow shock. The condition on the pulsar velocity for this tran- 
sition to occur while the rem nant is in the Sedov- Taylor phase is given by 
(Ivan der Swaluw et al.l . I2004J . and references therein) : 
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km s ^ , 



(15) 



a relation showing a strikingly weak dependence on the physical parameters. 
A significant fraction (30-4 0% depending on th e veloc ity distribution model) 
of the pulsars compiled by Arzoumanian et al.l (120021 ) satisfy this condition. 



van der Swaluw et al.l (120031 ) showed that once the pulsar reaches the edge of 
the remnant, its Mach number is M.psr — 3.1. Subsequently, the pulsar moves 
through the ISM where its velocity corresponds to a hypersonic Mach number 
typically on the order of 10^. 

The most f amous example of a PWN in this stage of evolution is the Gui- 
tar nebula (ICordes et al.l . Il993l . see Fig. 5), so named because of its cometary 
neck con necting to a nearly spherical b ubble. Numerous other examples are 
shown in iKargaltsev and Pavlovl (120081 ). A case of particular import to this 
work is that of the X-ray emission associated with PSR1929-I-10 (see Fig. 6). 
Wang et al.l (Il993l ) posited that the morphology is due to a relativistic back- 
fiow behind the pulsar, a suggestion that has gone unconfirmed for realistic 
wind Lorentz factors and pulsar velocities, and was a prime motivator for 
this project. The simulations in this section directly probe the morphology 
and interior structure of PWNe during this phase, motivate how the shape of 
the Guitar nebula persists, \ yithout resor t ing t o tailored ISM geometry, and 
confirm the interpretation of IWang et al.l (119931 ). 




» 1 




Fig. 5. A 1995 Hale Telescope Ha image of the Guitar Nebula (20 angstrom filter 
at 6564 angstro ms). The cometary neck conne cting to a spherical bubble are clearly 
evident. Credit: Chatterjee and Cordes ( 20021 ). 



5.2 Identifying suitable input parameters 



The outfiow streams relativistically into the ambient medium generating a 
strong shock. We derive a value for the outfiow pressure, Po, from the assump- 
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Fig. 6. -RO SAT X-ray s urface brightness in the field of PSR1929-I-10 showing the 
X-ray tail. Wang et al. ( 19931 ) suggested that the X-ray morphology is due to the 
acceleration of particles be hind the pulsar for ming a relativistic backflow. North is 



up and East is left. Credit: Wang et al. ( 19931 ) 



tion that the outflow is interacting with the ambient medium requiring that 
the momentum flux be comparable on either side of this shock; if the fluxes 
were not comparable, then either the ambient flow or outflow would dominate 
and the problem would be uninteresting. The momentum flux of the non- 
relativistic ambient medium and ultra-relativistic outflow are, respectively: 



FM,a = naVl +Pa, 

FM,o = loieo+Po)vl+Po- 

For an ultra-relativistic outflow, Po ^ no ^ Co ^ Spo, and Vo — > 1, and, for 
the ambient medium, UaV^ 3> Pa- Applying these conditions, and noting that 
7oPo > Po, gives: 



Po~n, l^-^j ~ 10-'' for 7„ = 10^ = 1. 

We are then free to pick any Uo meeting the conditions of a light, relativistic 
outflow, i.e. Ha, Po ^ Uo- This condition is motivated by the fact that the 
flow is very fast (7 > 10^). Well below the length scale of this study, the flow 
will be stabilized by the strong magnetic field, synchrotron cooling will be 
strong, and adiabatic losses due to expansion across the orders of magnitude 
in scale between the pulsar and the termination shock will sap internal energy. 
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This will conspire to effectively stop energy from being converted into thermal 
motions. Under those conditions, the flow might be cold. However, on the scale 
of the termination shock the field is much weaker meaning far less synchrotron 
cooling, instability is less inhibited, and interesting evolution will occur over 
fewer adiabatic-loss scale lengths. Indeed, the flow might well be influenced by 
waves generated both upstream and downstream of the shock(s). The scales 
relevant to t his work corr e spond to the region of "hot" post-termination shock 
plasma (e.g. iBucciantinil . |2008| ) within a PWN. Therefore, a hot flow seems 
significantly more plausible than does a cold flow and we select Uo = 10~'po, 
3 < / < 6. This clearly satisfies Po ^ and one may verify it satisfies Ua ^ rio 
by noting that the equation for po implies ^ po since 7^ ^ f ^ for the flows 
of interest here. 



5.3 A relativistic backflow 



Fig. 7 shows a simulation of a 70 = 10^ outflow interacting with an ambient 
flow with velocity Va = 0.00583(~1750 km s~^). The outflow pressure was cal- 
culated for an ambient-flow velocity of 500 km s~^ in order to match the typical 
value for pulsars in general. The outflow originates inside the circular region 
to the left of the evolving structure and the ambient flow streams in along the 
left edge of the computational domain. Fig. 8 shows the limited extent of the 
refined grid, supporting the choice of a maximum number of refinement level 
of Ljnax = 1- Recall the Ha image of the Guitar Nebula (see Fig. 5), a well- 
known pulsar wind nebula with the most rapidly movi ng pulsar ever observed , 



with a transverse velocity of (1.7±0.4) x 10'^ km s ^ (IChatterjee and Cordes 



2OO2I ). The simulation qualitatively resembles the nebula. This result consti- 
tutes compelling motivation for the conclusion that interstellar-medium flows 
set up by the space motion of pulsars can indeed produce "cometary" nebulae. 

We believe this simulation to be the first demonstrating asymmetry arising 
from a spherically-symmetric, light, ultra-relativistic flow interacting with a 
dense, slow ambient flow. The lines labeled "1" and "2" on the density map 
in Fig. 7 mark one-dimensional cuts (hereafter "cut-hl"and "cut-h2", respec- 
tively) made to probe the state of the simulation. Cut-hl spans the entire 
structure while cut-h2 spans the interior space occupied by the pressure en- 
hancements clearly visible in the pressure map. Fig. 9 shows the values of the 
flow parameters along these cuts. These plots clearly show the outer bounding 
Shockwave represented by the red boundary in the density map as well as a 
series of weaker internal on-axis shocks visible in the pressure map. The x- 
component of the flow velocity shows that a relativistic back flow harboring a 
se ries of weak shocks has arisen down stream. This validates the interpretation 



bv lWang et all (11993f ) of the origin of the X-ray trail behind PSR1929+10, and 



demonstrates the ability of the refined solver to elucidate the internal struc- 
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Fig. 7. An 871,200-iteration simulation of a light, ultra-relativistic outflow interact- 
ing with a dense, slow ambient flow. The input parameters are: Va = 0.00583 (~1750 
km s-i), M = 300, n„ = 1, 7^ = 10^ Po= 7xl0-^^ and = W'^Po. The upper 
and lower panels show a linear color map of the rest-frame pressure and lab-frame 
mass density, respectively. Both have been reflected along the symmetry axis. The 
outflow originates within the circular region to the left of the evolving structure and 
the ambient flow streams in along the left edge of the domain. The lines labelled 
"1" and "2" are 1-D data cuts (hereafter "cut-hl"and "cut-h2", respectively) with 
flow parameters plotted in Fig. 9. 
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Fig. 8. Plotted in red overlaying the pressure map for the simulation shown in Fig. 7 
is the refined grid at level L=l. The bottom half of the map is a refiection of the top 
half and has the same refined grid even though it is not shown. Note that the red 
lines trace the outlines of the meshes of refined cells, but not the cells themselves. 
While the boundary shock is well-refined, the axial shocks within the nebula are not 
refined at all. Flagging is determined by the largest difference in R between adjacent 
cells for all cells at level L. The refinement only follows the boundary shock because 
R differences inside the nebula are small compared to the difference between the 
nebula and the ambient medium. We will investigate refinement fiagging in more 
detail in a future study. 
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X (grid coordinates) 

ture of diffuse, ultra-relativistic pulsar wind nebulae which is often difficult to 
observe directly. 

It is noteworthy that the termination shock of the wind is not evident in the 
simulation discussed above. This is due to numerical shocking of the wind as 
it emerges from the on-axis hemisphere, as follows. Consider the cells depicted 
in Fig. 10. Let the angle of the line connecting the center of the hemisphere 
and the center of a cell 1, 2, 3, or 4 be 9i, i = 1,2, 3, 4. Since we have taken the 
pulsar wind to be spherically symmetric as it emerges from the hemesphere, 
we may calculate the relative flow velocities Avi2 and Af34 (normalized to the 
speed of light) at the centers of cells 1 & 2 and 3 & 4: 
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Fig. 9. The run of the laboratory-frame mass (R), momentum (M), and total energy 
(E) densities, rest-frame mass (n), and total energy (e) densities, and pressure (p), 
Lorentz factor (7), x- and y-components of the flow velocity (vx, Vy), the flow 
velocity (v), sound speed (cg), and generalized Mach number (Ai) along (a) cut-hl 
and (b) cut-h2 in Fig. 7. 
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Fig. 10. Schematic geometric cell layout pertaining to numerical shocking of the 
pulsar wind. The arc represents the on-axis hemisphere with radius 37.5 fine cells. 
Cell 1 is on-axis and is centered at 41.5 fine cells from the center of the hemisphere 
(relative center coordinates (x,y) = (41.5,0.5). The center coordinates of cells 2, 3, 
and 4 are (41.5,1.5), (29.5,29.5), and (29.5, 30.5) respectively. 

North 




East 



This shows that the relative velocity between vertically adjacent on-axis cells 
just outside the hemisphere is supersonic relative to the pulsar outflow sound 
speed of 0.57 (for the parameters relevant to Fig. 7). Thus, the wind near 
the axis shocks immediately and is thermalized producing a post-termination 
shock flow. Given that at early times the wind shows no deviation from spheri- 
cal symmetry (see Fig. 11), it is clear that this asymmetric numerical shocking 
of the wind is smeared out by the interaction with the ambient flow and does 
not impact the global evolution of the simulation. 

Additional refinement levels, perhaps needed only at early simulation times, 
will mitigate the numerical shocking issue. However, since tests have shown 
such shocking is present with 3 levels, and the significant results discussed 
below were possible with 2 (i.e. L^ax — 1), we leave explorations of addi- 
tional refinement levels to future studies. When these studies result in un- 
shocked, ultra-relativistic wind fiows entering the computational domain, and 
the resolution of the termination shock, we will perform new shock-tube tests. 
However, while the refined REST_FRAME routine is essential for proper han- 
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dling of the 7 >> 1 outflow, in the simulations presented here there are no 
structures involving Lorent z f actors higher than t hose previously explored by 
Duncan and Huehesl (Il994t ) fc lHuehes et~all tOO'j } with the RHLLE solver for 
which shock-tube validation was performed. All structures referenced below 
originated in the computational domain where the the tried-and-true Newton- 
Raphson iterative solver was toggled into action (recall §4). Therefore, we 
procede with flrm confldence rooted in the previous shock-tube tests. 



Fig. 11. A time sequence of linear pressure maps for the simulation shown in Fig. 7. 
The sequence indicates that the appearance of the inflection point is preceded by a 
pressure drop inside the nebula. Note the finer time steps between 240K and 360K 
iterations and that the color map is relative to the minimum and maximum for each 
plot individually. However, the minimum is the same and the maximum is similar 
for all plots, so the variation is minimal. 



(a) 10,000 iterations. 



(b) 120,000 iterations. 



(c) 240,000 iterations. 





(d) 300,000 iterations. (e) 320,000 iterations. 



(f) 360,000 iterations. 




(g) 480,000 iterations. (h) 600,000 iterations. 



(i) 720,000 iterations. 



6 Discussion 



The physics behind the formation of the structure observed in Fig. 7 is as 
follows. The wind streams outward and sweeps up ambient material which 
drives pressure waves (weakly at flrst) into the shocked wind. As the nebula 
expands, the pressure inside decreases. Fig. 11 shows a series of pressure maps 
comprising a time development sequence for the simulation depicted in Fig. 7. 
The sequence clearly indicates that once the pulsar crosses the boundary of 
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the initially spherical nebula, an inflection point develops along the leading 
edge at approximately 45° from the measured from W to nE]. The 

pressure waves are intensifled and propagate to the axis and reflect, leading 
to the formation of a relativistic backflow harboring internal Shockwaves rem- 
iniscent of shock diamonds. The fact that the backflow does not develop until 
after the inflection point supports this picture. The internal Shockwaves, in 
turn, thermalize energy, allowing the flow to expand and inflate the trailing 
spherical bubble. As the bubble inflates, it "pinches" the inflection point en- 
hancing the cuspy shape, maintaining the pressure-wave influx that sets up 
the energy-thermalizing backflow responsible for inflating the bubble. Such a 
feedback cycle is relevent to the Guitar nebula even though the pulsar was 
not born at the center of the trailing bubble - given its proper motion, the 
pul sar moves a di s tance corresponding to the entire nebula in less than 500 



yr (IRomani et al.l . Il997l ). a time orders of magnitude too short for the age 
of a pulsar powering a bow-shock nebula - because it explains how the bub- 
ble persists. Such a scenario is analogous to the formation of structure in 
relativistic galactic jets, where the evolution is driven by Kelvin-Helmholtz 
modes long the contact surface th at separates the sho cked ambient medium 



from the shocked jet material (e.g. iHughes et al.l . |2002| ) 



The evolution of the Guitar-like shape is rather sensitive to the choice of pa- 
rameters. As Tab. 3 shows, the appearance of the inflection point marking the 
onset of the formation of the "neck" of the Guitar takes a signiflcantly larger 
number of computational iterations as the ambient-flow velocity decreases. 
This is expected as the asymmetry of the nebula should evolve more slowly in 
this scenario: a decreased ambient-flow velocity is equivalent to a lower pulsar 
space motion and so it takes more time for the pulsar to reach the nebula 
boundary thereby delaying the formation of the inflection point. If the inflec- 
tion point is not induced while the expanding nebula is small enough such 
that the ensuing neck is of signiflcant scale, then no Guitar-like structure will 
be apparent. If a pulsar velocity of 1500 km s~^, 1250 km s~^, or 1000 km s~^ 
is required for observable Guitar-like mo rphology to arise, then the velocity 
distribution of lArzoumanian et al.l (|2002| ) implies that <5%, 7-8%, or ~15% 



of radio pulsars, respectively, have the possibility of developing such features 
depending on the nature of their ambient environment. 



7 Conclusion 



We discussed the application of an existing special relativistic, hydrodynamic 
(SRHD) primitive- variable recovery algorithm to ultra-relativistic flows (Lorentz 

^ This is sensible as it is the location where the wind velocity transitions from 
having its largest component at 180° to the inflow direction to having it at 90°. 
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Table 3 

The dependence of the Guitar-hke inflection point on the number of simulation iter- 
ations. As expected, the higher the ambient-flow velocity, the sooner the inflection 
point d('V(^lops due lo tlic iucixviscd rale al which aml)i('iil malcrird is s\vcp(-up. 



Wind Lorentz factor 


ambient-flow velocity 


Iterations until inflection 


(unitless) 


(km s^^) 


(10^) 










1750 


30 




1500 


39 




1250 


54 




1000 


81 




750 


1 m Gfif^n a + 7 A. 


10^ 








5500 


4 




4250 


6 




3000 


12 




1750 


30 



factor, 7, of 10^-10®) and the refinement necessary for the numerical velocity 
root finder to work in this domain. We found that the velocity quartic, Q{v), 
exhibits dual roots in the physical velocity range that move progressively closer 
together for larger 7 leading to a divide by zero and the failure of the Ncwton- 
Raphson iteration method employed by the existing primitives algorithm. Our 
solution was to recast the quartic to be a function, Q{j), of 7. We demon- 
strated that (5(7) exhibits only one physical root. However, Newton-Raphson 
iteration also failed in this case at high 7, due to the extreme slope of the 
quartic near the root, necessitating the use there of an analytical numerical 
root finder. 

Our timing analysis indicated that using Q{'y) with the 8-byte analytical root 
finder increased run time by only 24% compared to using Q{v) with the 8-byte 
iterative root finder (based on 10 trial runs), while using Q{j) with the 16-byte 
analytical root finder ballooned run time by a factor of approximately 400. The 
iterative root finder is accurate to order 10~^ for a sizable region of parameter 
space corresponding to Lorentz factors on the order of 10^ and smaller. There- 
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fore, we implemented a computational switch that checks the values of M/E 
and R/ E and calls the iterative or analytical root finder accordingly thereby 
creating a hybrid primitives recovery algorithm called REST_FRAME. 

In addition, our exploration of parameter space suggests that the discriminant 
of the cubic resolvent (as defined by Eqn. 9 in §3.1) will always be positive for 
physical fiows. Therefore, we did not include code for negative discriminants 
in our routine. Formal proof remains elusive, however, leaving potential for 
future work. 

We have shown that REST_FRAME is capable of calculating the primitive 
variables from the conserved variables to an accuracy of at least 0(10""^) 
for Lorentz factors up to 10^ with significantly better accuracy for Lorentz 
factors < 10^, and shghtly worse (order 10~^) for a small portion of the space 
corresponding to the highest Lorentz factors. We traced the degradation in 
accuracy for larger Lorentz factors to the effect of subtractive cancellation. 
Past studies have shown that an accuracy of order 10~^ is capable of robustly 
capturing hydrodynamic structures. We have applied the refined solver to an 
ultra-relativistic problem and have shown that it is capable of reproducing 
observed structures and is well-suited to our study of the internal structure of 
diffuse pulsar wind nebulae. 

Our main conclusions are as follows: 

• Relativistic, hydrodynamic simulations have shown that the relatively slow, 
dense ISM fiow resulting from the space motion of a pulsar can set up an 
interaction with the extremely light, ultra-relativistic pulsar wind leading to 
an asymmetric nebula with a morphology reminiscent of the Guitar nebula. 

• Simulations have validated the interpretation that a relativistic backfiow 
behind PSR1929+10 is responsible for the X-ray morphology. Results fur- 
ther show that the backfiow can harbor a series of internal Shockwaves that 
infiates a nebular bubble, and that the bubble provides positive feedback to 
the backfiow, explaining how the Guitar bubble persists. 

• The evolution of the bubble/backfiow structure is sensitive to the choice 
of input parameters justifying a future series of simulation runs that will 
determine what pulsar velocities and wind/ISM density ratios are required 
for the bubble/backfiow feedback loop to arise. 
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